Fractional Diffusion based on Riemann-Liouville Fractional Derivatives

A fractional diffusion equation based on Riemann-Liouville fractional derivatives is solved exactly. The initial values are given as fractional integrals. The solution is obtained in terms of H-functions. It differs from the known solution of fractional diffusion equations based on fractional integrals. The solution of fractional diffusion based on a Riemann-Liouville fractional time derivative does not admit a probabilistic interpretation in contrast with fractional diffusion based on fractional integrals. While the fractional initial value problem is well defined and the solution finite at all times its values for t → 0 are divergent. J.Phys.Chem B, vol. 104, page 3914, (2000) Anomalous subdiffusive transport appears to be a universal experimental phenomenon [1–3]. Examples occur in widely different systems ranging from amorphous semiconductors [4,5] through polymers [6–8] and composite heterogeneous films [9] to porous media [10,11]. Theoretical investigations into anomalous diffusion and continuous time random walks have been a major focus of H. Schers research for many years [12,13,4,10]. The purpose of this paper is to discuss a theoretical approach based on the replacement of the time derivative in the diffusion equation with a derivative of noninteger order (fractional derivative). Many investigators have proposed the use of fractional time derivatives for subdiffusive transport on a purely mathematical or heuristic basis [14–19,8,20]. From the perspective of theoretical physics this proposal touches upon fundamental principles such as locality, irreversibility and invariance under time translations because fractional derivatives are nonlocal operators that are not invariant under time reversal [21]. These issues are generally avoided in heuristic and mathematical proposals, but were discussed recently in the context of long time limits and coarse graining [21]. It was found that fractional time derivatives with orders between 0 and 1 may generally appear as infinitesimal generators of a coarse grained macroscopic time evolution [20–24]. Differential equations involving fractional derivatives raise a second basic problem, related to the first, that will be the focus of this paper. The second problem is whether to replace the integer order derivative by a Riemann-Liouville, by a Weyl, by a Riesz, by a Grünwald or by a Marchaud fractional derivative (see [25–28] for definitions of these different derivatives). Different authors have introduced different derivatives depending on the physical situation [21,29–31]. Given the basic objective of introducing fractional derivatives into the diffusion equation the present paper will be concerned with the equation Dα0+f(r, t) = Cα∆f(r, t) (1) where f(r, t) denotes the unknown field and Cα denotes the fractional diffusion constant with dimensions [cm/s]. The fractional derivative operator, denoted as Dα0+, is the Riemann1 Liouville derivative of order α and with lower limit a ∈ R. It is defined as [25,28] (Dαa+f)(x) = d dx (I a+ f)(x) (2) where (Ia+f)(x) = 1 Γ(α) ∫ x a (x− y)f(y) dy (3) is the Riemann-Liouville fractional integral with order α and lower limit a. Although many authors have investigated fractional diffusion problems [14,15,32,16,19,33] it seems that equation (1) has not been solved previously. In fact it was recently questioned whether an approach using eq. (1) is consistent [34,29]. It is the purpose of this paper to solve eq. (1) exactly thereby establishing its consistency for appropriate initial conditions. Let me emphasize that eq. (1) differs from the popular equation introduced and solved in [16]. The latter equation is obtained by first rewriting the diffusion equation in integral form as f(r, t) = f0δ(r) + C1 ∫ t 0 ∆f(r, t) dt (4) where C1 is the usual diffusion constant, δ(r) is the Dirac measure at the origin, and where the initial condition f(r, 0) = f0δ(r) has been incorporated. Then the integral on the right hand side is replaced by a fractional Riemann-Liouville integral to arrive at the fractional integral form f(r, t) = f0δ(r) + Cα Γ(α) ∫ t 0 (t− t)∆f(r, t) dt = f0δ(r) + Cα(I0+∆f)(r, t) (5a) or, upon differentiating both sides, at ∂ ∂t f(r, t) = Cα(D 1−α 0+ ∆f)(r, t) (5b) where Cα is again a fractional diffusion constant. For α = 1 this reduces to eq. (4). The exact solution of eq. (5) is known and given by eq. (22) below. In [35,36] it was shown that eqs. (5) have a rigorous relation with continuous time random walks of the kind investigated frequently by Harvey Scher [12,13,4,10]. More precisely, eq. 2 (1) was found to correspond exactly to a continuous time random walk with the long tailed waiting time density ψ(t;α, τ0) = 1 τ0 ( t τ0 α−1 Eα,α ( − t α τ 0 ) (6) where τ0 is a time constant. Here Ea,b(x) denotes the generalized Mittag-Leffler function defined by Ea,b(x) = ∞ ∑ k=0 x Γ(ak + b) (7) for all a > 0 and b ∈ C. For α = 1 this reduces to an exponential waiting time density. For 0 < α < 1 these waiting time densities have a long tail decaying as ψ(t) ∼ t for t → ∞. Interestingly ψ(t) ∼ t diverges algebraically for t → 0. It follows from refs. [22,23,20,24,21] that among the the waiting time densities with long tails the densities ψ(t;α, τ0) represent important universality classes for continuous time random walks. Note that eqs. (5) and (1) are not equivalent. The difference between eqs. (5) and (1) has to do with the initial conditions. An appropriate inital condition is found by analysing the stationary case. One finds that the fractional integral I (1−α) 0+ f(r, 0+) = f0,αδ(r) (8) is preserved during the time evolution. This is a nonlocal initial condition. It implies the divergence of f(r, t) as t→ 0, as is characterstic for fractional stationarity [20,24]. Equation (1) with initial condition (8) can be solved exactly by Fourier-Laplace techniques. Let the Fourier transformation be defined as F {f(r)} (q) = ∫ R ef(r)dr. (9) Fourier and Laplace transformation of eq. (1) now yield f(q, u) = f0,α Cαq + uα . (10) Inverting the Laplace transform gives f(q, t) = f0,α t Eα,α(−Cαqt). (11) 3 Setting q = 0 shows that f(r, t) cannot be a probability density because its normalization would depend on t. Hence eq. (1) does not admit a probabilistic interpretation contrary to eq. (5). To obtain f(r, t) it is advantageous to first invert the Fourier transform in eq. (10) and only later the Laplace transform. The Fourier transform may be inverted by noting the formula [37] (2π) ∫ e ·r ( |r| m 1−(d/2) K(d−2)/2 (m|r|) dr = 1 q2 +m2 (12) which leads to f(r, u) = f0,α(2πCα) −d/2 ( r √ Cα )1−(d/2) uK(d−2)/2 ( ru √ Cα ) (13) with r = |r|. To invert the Laplace transform it is convenient to use the relation M{f(t)} (s) = M{L{f(t)} (u)} (1− s) Γ(1− s) (14) between the Laplace transform and the Mellin transform M{f(t)} (s) = ∫ ∞ 0 tf(t) dt. (15) of a function f(t). Setting A = r/ √ Cα, λ = α/2, ν = (d − 2)/2 and μ = α(d − 2)/4 and using the general relation M{xg(bx)} (s) = 1 p bg ( s+ q p ) (b, p > 0) (16) leads to M{f(r, u)} (s) = f0,α λ (2πCα) AAM{Kν(u)} ((s+ μ)/λ) . (17) The Mellin transform of the Bessel function reads [38] M{Kν(x)} (s) = 2Γ ( s+ ν 2 )